Calculate Schoener’s overlap and Levin’s diversity index and fit brms models
# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(bayesplot)
library(tidybayes)
library(brms)
# To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/cod_flounder_diets_spatial_cache/html")
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf()
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
theme_clean <- function() {
theme_minimal(base_family = "Barlow Semi Condensed") +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
strip.text = element_text(family = "BarlowSemiCondensed-Bold",
size = rel(1), hjust = 0),
strip.background = element_rect(fill = "grey80", color = NA))
}
# This is just to add the density covariates to the schoener overlap models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> time_period = col_character(),
#> quarter_fact = col_character(),
#> pred_weight_source = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
cod <- cod %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(year > 2014) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 102 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 1,666 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 7,108 unique values and 0% NA
#> new variable 'depth_sc' (double) with 102 unique values and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed
#> drop_na: no rows removed
# Now read stomach data (1 row 1 stomach)
cod_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_full_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "COD") %>% filter(year > 2014) %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 4,936 rows (58%), 3,593 rows remaining
#> filter: no rows removed
fle_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_full_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "FLE") %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: no rows removed
# Plot data in space
plot_map_raster +
geom_point(data = fle_prey, aes(x = X * 1000, y = Y * 1000, color = "FLE"), alpha = 0.5) +
geom_point(data = cod_prey, aes(x = X * 1000, y = Y * 1000, color = "COD"), alpha = 0.5) +
theme_plot() +
facet_wrap(~sub_div)
cod_prey <- cod_prey %>% filter(lat < 58)
#> filter: no rows removed
fle_prey <- fle_prey %>% filter(lat < 58)
#> filter: no rows removed
# I will also make a new area category, pooling 24 and 25, and 27-8, making it a southwest (coastal), north (coastal) and offshore (26). Note this is only for the cod_prey and fle_prey data, which we use for indices. The other stomach analysis is spatial, but here we need enough samples to calculate indices properly
cod_prey <- cod_prey %>% mutate(area = ifelse(sub_div %in% c(24, 25), "24-25", NA),
area = ifelse(sub_div %in% c(27, 28), "27, 28", area),
area = ifelse(sub_div == 26, "26", area))
#> mutate: new variable 'area' (character) with 3 unique values and 0% NA
fle_prey <- fle_prey %>% mutate(area = ifelse(sub_div %in% c(24, 25), "24-25", NA),
area = ifelse(sub_div %in% c(27, 28), "27, 28", area),
area = ifelse(sub_div == 26, "26", area))
#> mutate: new variable 'area' (character) with 3 unique values and 0% NA
# Reformat data
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, area) %>%
summarise(fle_stomach_content = sum(value)) %>%
arrange(name, year, area, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, area, sep = "_")) %>%
group_by(id) %>%
mutate(fle_stomach_content_tot = sum(fle_stomach_content),
fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x32, now 34880x18]
#> group_by: 4 grouping variables (name, year, quarter, area)
#> summarise: now 352 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 22 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 22 unique values and 0% NA
#> new variable 'fle_stomach_content_prop' (double) with 199 unique values and 0% NA
#> ungroup: no grouping variables
# This should amount to the number of unique prey
number_of_prey <- fle_prey_long %>% group_by(id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 351 rows (>99%), one row remaining
number_of_prey
#> # A tibble: 1 × 1
#> n
#> <int>
#> 1 16
test_id <- head(fle_prey_long$id, 1)
fle_prey_long %>% filter(id == test_id) %>% as.data.frame()
#> filter: removed 336 rows (95%), 16 rows remaining
#> name year quarter area fle_stomach_content id
#> 1 amphipoda_tot 2015 4 24-25 7.64 2015_4_24-25
#> 2 bivalvia_tot 2015 4 24-25 291.27 2015_4_24-25
#> 3 clupea_harengus_tot 2015 4 24-25 0.00 2015_4_24-25
#> 4 clupeidae_tot 2015 4 24-25 0.00 2015_4_24-25
#> 5 gadiformes_tot 2015 4 24-25 0.00 2015_4_24-25
#> 6 gadus_morhua_tot 2015 4 24-25 0.00 2015_4_24-25
#> 7 gobiidae_tot 2015 4 24-25 0.00 2015_4_24-25
#> 8 mysidae_tot 2015 4 24-25 0.00 2015_4_24-25
#> 9 non_bio_tot 2015 4 24-25 6.59 2015_4_24-25
#> 10 other_crustacea_tot 2015 4 24-25 4.05 2015_4_24-25
#> 11 other_pisces_tot 2015 4 24-25 0.00 2015_4_24-25
#> 12 other_tot 2015 4 24-25 3.44 2015_4_24-25
#> 13 platichthys_flesus_tot 2015 4 24-25 0.00 2015_4_24-25
#> 14 polychaeta_tot 2015 4 24-25 2.82 2015_4_24-25
#> 15 saduria_entomon_tot 2015 4 24-25 45.85 2015_4_24-25
#> 16 sprattus_sprattus_tot 2015 4 24-25 0.00 2015_4_24-25
#> fle_stomach_content_tot fle_stomach_content_prop
#> 1 361.66 0.021124813
#> 2 361.66 0.805369684
#> 3 361.66 0.000000000
#> 4 361.66 0.000000000
#> 5 361.66 0.000000000
#> 6 361.66 0.000000000
#> 7 361.66 0.000000000
#> 8 361.66 0.000000000
#> 9 361.66 0.018221534
#> 10 361.66 0.011198363
#> 11 361.66 0.000000000
#> 12 361.66 0.009511696
#> 13 361.66 0.000000000
#> 14 361.66 0.007797379
#> 15 361.66 0.126776530
#> 16 361.66 0.000000000
# Now cod
# colnames(cod_prey)
cod_prey_long <- cod_prey %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, area) %>%
summarise(cod_stomach_content = sum(value)) %>%
arrange(name, year, area, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, area, sep = "_")) %>%
group_by(id) %>%
mutate(cod_stomach_content_tot = sum(cod_stomach_content),
cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x32, now 57488x18]
#> group_by: 4 grouping variables (name, year, quarter, area)
#> summarise: now 416 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 26 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 26 unique values and 0% NA
#> new variable 'cod_stomach_content_prop' (double) with 299 unique values and 0% NA
#> ungroup: no grouping variables
unique(is.na(fle_prey_long))
#> name year quarter area fle_stomach_content id
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> fle_stomach_content_tot fle_stomach_content_prop
#> [1,] FALSE FALSE
unique(is.na(cod_prey_long))
#> name year quarter area cod_stomach_content id
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> cod_stomach_content_tot cod_stomach_content_prop
#> [1,] FALSE FALSE
fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, -9)) %>% filter(fle_stomach_content_prop == -9)
#> mutate: no changes
#> filter: removed all rows (100%)
#> # A tibble: 0 × 8
#> # … with 8 variables: name <chr>, year <dbl>, quarter <dbl>, area <chr>,
#> # fle_stomach_content <dbl>, id <chr>, fle_stomach_content_tot <dbl>,
#> # fle_stomach_content_prop <dbl>
fle_prey_long <- fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, 0))
#> mutate: no changes
# Calculate Levin niche index
levin <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_stomach_content_prop) %>%
drop_na(cod_stomach_content_prop) %>%
group_by(year, quarter, area) %>%
summarise(levin_cod = (1/(number_of_prey$n-1)) * (((1/sum(cod_stomach_content_prop^2))) - 1),
levin_fle = (1/(number_of_prey$n-1)) * (((1/sum(fle_stomach_content_prop^2))) - 1)) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "area", "id")
#> left_join: added 3 columns (fle_stomach_content, fle_stomach_content_tot, fle_stomach_content_prop)
#> > rows only in x 64
#> > rows only in y ( 0)
#> > matched rows 352
#> > =====
#> > rows total 416
#> drop_na: no rows removed
#> drop_na: removed 64 rows (15%), 352 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, area)
#> summarise: now 22 rows and 5 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
levin
#> # A tibble: 22 × 5
#> year quarter area levin_cod levin_fle
#> <dbl> <dbl> <chr> <dbl> <dbl>
#> 1 2015 4 24-25 0.216 0.0335
#> 2 2015 4 26 0.111 0.152
#> 3 2015 4 27, 28 0.232 0.0623
#> 4 2016 1 24-25 0.154 0.174
#> 5 2016 1 26 0.165 0.0352
#> 6 2016 1 27, 28 0.0736 0.112
#> 7 2016 4 24-25 0.112 0.0231
#> 8 2016 4 26 0.137 0.0292
#> 9 2016 4 27, 28 0.0780 0.122
#> 10 2017 1 24-25 0.146 0.0602
#> # … with 12 more rows
# Quickly check data to determine which distribution to use
levin %>%
ungroup() %>%
count(levin_cod == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now one row and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with one unique value and 0% NA
#> # A tibble: 1 × 3
#> `levin_cod == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 22 1
levin %>%
ungroup() %>%
count(levin_fle == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now one row and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with one unique value and 0% NA
#> # A tibble: 1 × 3
#> `levin_fle == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 22 1
# Fit beta models, so few zeroes and no 1's
levin %>% arrange(desc(levin_cod)) %>% dplyr::select(levin_cod)
#> # A tibble: 22 × 1
#> levin_cod
#> <dbl>
#> 1 0.246
#> 2 0.232
#> 3 0.216
#> 4 0.188
#> 5 0.165
#> 6 0.159
#> 7 0.154
#> 8 0.149
#> 9 0.146
#> 10 0.137
#> # … with 12 more rows
levin %>% arrange(desc(levin_fle)) %>% dplyr::select(levin_fle)
#> # A tibble: 22 × 1
#> levin_fle
#> <dbl>
#> 1 0.184
#> 2 0.174
#> 3 0.166
#> 4 0.153
#> 5 0.152
#> 6 0.129
#> 7 0.122
#> 8 0.114
#> 9 0.112
#> 10 0.0993
#> # … with 12 more rows
# Final polish of data before feeding into models (species, not size-based indicies)
levin <- levin %>%
mutate(levin_cod2 = ifelse(levin_cod == 0, 0.0001, levin_cod),
levin_fle2 = ifelse(levin_fle == 0, 0.0001, levin_fle),
year_f = as.factor(year),
quarter_f = as.factor(quarter),
area_f = as.factor(area))
#> mutate: new variable 'levin_cod2' (double) with 22 unique values and 0% NA
#> new variable 'levin_fle2' (double) with 22 unique values and 0% NA
#> new variable 'year_f' (factor) with 6 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> new variable 'area_f' (factor) with 3 unique values and 0% NA
# Reformat data to calculate Schoeners overlap index
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, ices_rect) %>%
summarise(fle_stomach_content = sum(value)) %>%
arrange(name, year, ices_rect, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>%
group_by(id) %>%
mutate(fle_stomach_content_tot = sum(fle_stomach_content),
fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x32, now 34880x18]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,296 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 81 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 81 unique values and 0% NA
#> new variable 'fle_stomach_content_prop' (double) with 445 unique values and 1% NA
#> ungroup: no grouping variables
# This should amount to the number of unique prey
number_of_prey <- fle_prey_long %>% group_by(id) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 1,295 rows (>99%), one row remaining
number_of_prey
#> # A tibble: 1 × 1
#> n
#> <int>
#> 1 16
test_id <- head(fle_prey_long$id, 1)
fle_prey_long %>% filter(id == test_id) %>% as.data.frame()
#> filter: removed 1,280 rows (99%), 16 rows remaining
#> name year quarter ices_rect fle_stomach_content
#> 1 amphipoda_tot 2015 4 40G4 1.30
#> 2 bivalvia_tot 2015 4 40G4 240.98
#> 3 clupea_harengus_tot 2015 4 40G4 0.00
#> 4 clupeidae_tot 2015 4 40G4 0.00
#> 5 gadiformes_tot 2015 4 40G4 0.00
#> 6 gadus_morhua_tot 2015 4 40G4 0.00
#> 7 gobiidae_tot 2015 4 40G4 0.00
#> 8 mysidae_tot 2015 4 40G4 0.00
#> 9 non_bio_tot 2015 4 40G4 5.78
#> 10 other_crustacea_tot 2015 4 40G4 2.12
#> 11 other_pisces_tot 2015 4 40G4 0.00
#> 12 other_tot 2015 4 40G4 2.27
#> 13 platichthys_flesus_tot 2015 4 40G4 0.00
#> 14 polychaeta_tot 2015 4 40G4 2.82
#> 15 saduria_entomon_tot 2015 4 40G4 0.07
#> 16 sprattus_sprattus_tot 2015 4 40G4 0.00
#> id fle_stomach_content_tot fle_stomach_content_prop
#> 1 2015_4_40G4 255.34 0.0050912509
#> 2 2015_4_40G4 255.34 0.9437612595
#> 3 2015_4_40G4 255.34 0.0000000000
#> 4 2015_4_40G4 255.34 0.0000000000
#> 5 2015_4_40G4 255.34 0.0000000000
#> 6 2015_4_40G4 255.34 0.0000000000
#> 7 2015_4_40G4 255.34 0.0000000000
#> 8 2015_4_40G4 255.34 0.0000000000
#> 9 2015_4_40G4 255.34 0.0226364847
#> 10 2015_4_40G4 255.34 0.0083026553
#> 11 2015_4_40G4 255.34 0.0000000000
#> 12 2015_4_40G4 255.34 0.0088901073
#> 13 2015_4_40G4 255.34 0.0000000000
#> 14 2015_4_40G4 255.34 0.0110440981
#> 15 2015_4_40G4 255.34 0.0002741443
#> 16 2015_4_40G4 255.34 0.0000000000
# Now cod
# colnames(cod_prey)
cod_prey_long <- cod_prey %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, ices_rect) %>%
summarise(cod_stomach_content = sum(value)) %>%
arrange(name, year, ices_rect, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>%
group_by(id) %>%
mutate(cod_stomach_content_tot = sum(cod_stomach_content),
cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x32, now 57488x18]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,696 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 106 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 106 unique values and 0% NA
#> new variable 'cod_stomach_content_prop' (double) with 784 unique values and 0% NA
#> ungroup: no grouping variables
unique(is.na(fle_prey_long))
#> name year quarter ices_rect fle_stomach_content id
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> [2,] FALSE FALSE FALSE FALSE FALSE FALSE
#> fle_stomach_content_tot fle_stomach_content_prop
#> [1,] FALSE FALSE
#> [2,] FALSE TRUE
unique(is.na(cod_prey_long))
#> name year quarter ices_rect cod_stomach_content id
#> [1,] FALSE FALSE FALSE FALSE FALSE FALSE
#> cod_stomach_content_tot cod_stomach_content_prop
#> [1,] FALSE FALSE
fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, -9)) %>% filter(fle_stomach_content_prop == -9)
#> mutate: changed 16 values (1%) of 'fle_stomach_content_prop' (16 fewer NA)
#> filter: removed 1,280 rows (99%), 16 rows remaining
#> # A tibble: 16 × 8
#> name year quarter ices_rect fle_stomach_cont… id fle_stomach_cont…
#> <chr> <dbl> <dbl> <chr> <dbl> <chr> <dbl>
#> 1 amphipoda… 2016 1 40G8 0 2016_… 0
#> 2 bivalvia_… 2016 1 40G8 0 2016_… 0
#> 3 clupea_ha… 2016 1 40G8 0 2016_… 0
#> 4 clupeidae… 2016 1 40G8 0 2016_… 0
#> 5 gadiforme… 2016 1 40G8 0 2016_… 0
#> 6 gadus_mor… 2016 1 40G8 0 2016_… 0
#> 7 gobiidae_… 2016 1 40G8 0 2016_… 0
#> 8 mysidae_t… 2016 1 40G8 0 2016_… 0
#> 9 non_bio_t… 2016 1 40G8 0 2016_… 0
#> 10 other_cru… 2016 1 40G8 0 2016_… 0
#> 11 other_pis… 2016 1 40G8 0 2016_… 0
#> 12 other_tot 2016 1 40G8 0 2016_… 0
#> 13 platichth… 2016 1 40G8 0 2016_… 0
#> 14 polychaet… 2016 1 40G8 0 2016_… 0
#> 15 saduria_e… 2016 1 40G8 0 2016_… 0
#> 16 sprattus_… 2016 1 40G8 0 2016_… 0
#> # … with 1 more variable: fle_stomach_content_prop <dbl>
fle_prey_long <- fle_prey_long %>% mutate(fle_stomach_content_prop = replace_na(fle_stomach_content_prop, 0))
#> mutate: changed 16 values (1%) of 'fle_stomach_content_prop' (16 fewer NA)
# Calculate Schoener index
schoener <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_stomach_content_prop) %>%
drop_na(cod_stomach_content_prop) %>%
group_by(year, quarter, ices_rect) %>%
summarise(schoener = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop)))) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect", "id")
#> left_join: added 3 columns (fle_stomach_content, fle_stomach_content_tot, fle_stomach_content_prop)
#> > rows only in x 432
#> > rows only in y ( 32)
#> > matched rows 1,264
#> > =======
#> > rows total 1,696
#> drop_na: no rows removed
#> drop_na: removed 432 rows (25%), 1,264 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 79 rows and 4 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
# Summarise cod and flounder data by ices_rect then add to diet data
colnames(cod)
#> [1] "FR_tot" "FR_sad"
#> [3] "FR_spr" "FR_her"
#> [5] "saduria_entomon_tot" "tot_prey_biom"
#> [7] "common_tot" "unique_pred_id"
#> [9] "year" "quarter"
#> [11] "time_period" "quarter_fact"
#> [13] "pred_weight_g" "pred_weight_source"
#> [15] "pred_cm" "predator_spec"
#> [17] "predcod_density" "predfle_density"
#> [19] "predcod_density_sc" "predfle_density_sc"
#> [21] "depth" "X"
#> [23] "Y" "lat"
#> [25] "long" "ices_rect"
#> [27] "sub_div" "cruise"
#> [29] "depth2_sc" "saduria_entomon_per_mass"
#> [31] "tot_prey_biom_per_mass" "depth_sc"
cod$year_rect_id <- paste(cod$year, cod$quarter, cod$ices_rect, sep = "_")
dens_sum <- cod %>%
drop_na(predfle_density) %>%
drop_na(predcod_density) %>%
group_by(year_rect_id) %>%
summarise(predfle_density_tot = sum(predfle_density),
predcod_density_tot = sum(predcod_density)) %>%
ungroup() %>%
mutate(predfle_density_tot_sc = (predfle_density_tot - mean(predfle_density_tot)) / sd(predfle_density_tot),
predcod_density_tot_sc = (predcod_density_tot - mean(predcod_density_tot)) / sd(predcod_density_tot))
#> drop_na: no rows removed
#> drop_na: no rows removed
#> group_by: one grouping variable (year_rect_id)
#> summarise: now 106 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'predfle_density_tot_sc' (double) with 106 unique values and 0% NA
#> new variable 'predcod_density_tot_sc' (double) with 106 unique values and 0% NA
schoener$year_rect_id <- paste(schoener$year, schoener$quarter, schoener$ices_rect, sep = "_")
schoener <- left_join(schoener, dens_sum)
#> Joining, by = "year_rect_id"
#> left_join: added 4 columns (predfle_density_tot, predcod_density_tot, predfle_density_tot_sc, predcod_density_tot_sc)
#> > rows only in x 0
#> > rows only in y (27)
#> > matched rows 79
#> > ====
#> > rows total 79
# Quickly check data to determine which distribution to use
schoener %>%
ungroup() %>%
count(schoener == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `schoener == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 76 0.962
#> 2 TRUE 3 0.0380
# Fit beta models, so few zeroes and no 1's
schoener %>% arrange(desc(schoener)) %>% dplyr::select(schoener)
#> # A tibble: 79 × 1
#> schoener
#> <dbl>
#> 1 0.624
#> 2 0.614
#> 3 0.565
#> 4 0.556
#> 5 0.553
#> 6 0.5
#> 7 0.349
#> 8 0.308
#> 9 0.222
#> 10 0.208
#> # … with 69 more rows
# Final polish of data before feeding into models (species, not size-based indicies)
schoener <- schoener %>%
mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
year_f = as.factor(year),
quarter_f = as.factor(quarter),
ices_rect_f = as.factor(ices_rect))
#> mutate: new variable 'schoener2' (double) with 77 unique values and 0% NA
#> new variable 'year_f' (factor) with 6 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA
# colnames(cod_prey)
cod_prey_long2 <- cod_prey %>%
mutate(size_group = ifelse(pred_cm > 30, "Large", "Small")) %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, ices_rect, size_group) %>%
summarise(cod_stomach_content = sum(value)) %>%
arrange(name, year, ices_rect, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, ices_rect, sep = "_")) %>%
group_by(id, size_group) %>% # Add size-group here
mutate(cod_stomach_content_tot = sum(cod_stomach_content),
cod_stomach_content_prop = cod_stomach_content
/ cod_stomach_content_tot) %>%
ungroup()
#> mutate: new variable 'size_group' (character) with 2 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> group_by: 5 grouping variables (name, year, quarter, ices_rect, size_group)
#> summarise: now 3,200 rows and 6 columns, 4 group variables remaining (name, year, quarter, ices_rect)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 106 unique values and 0% NA
#> group_by: 2 grouping variables (id, size_group)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 193 unique values and 0% NA
#> new variable 'cod_stomach_content_prop' (double) with 1,079 unique values and <1% NA
#> ungroup: no grouping variables
# This should amount to the number of unique prey
cod_prey_long2 %>% group_by(id, size_group) %>% mutate(n = n()) %>% ungroup() %>% distinct(n)
#> group_by: 2 grouping variables (id, size_group)
#> mutate (grouped): new variable 'n' (integer) with one unique value and 0% NA
#> ungroup: no grouping variables
#> distinct: removed 3,199 rows (>99%), one row remaining
#> # A tibble: 1 × 1
#> n
#> <int>
#> 1 16
# Split by size-group
cod_prey_long_small <- cod_prey_long2 %>% filter(size_group == "Small") %>%
rename("cod_stomach_content_prop_small" = "cod_stomach_content_prop") %>%
dplyr::select(name, id, cod_stomach_content_prop_small)
#> filter: removed 1,632 rows (51%), 1,568 rows remaining
#> rename: renamed one variable (cod_stomach_content_prop_small)
cod_prey_long_large <- cod_prey_long2 %>% filter(size_group == "Large") %>%
rename("cod_stomach_content_prop_large" = "cod_stomach_content_prop") %>%
dplyr::select(name, id, cod_stomach_content_prop_large)
#> filter: removed 1,568 rows (49%), 1,632 rows remaining
#> rename: renamed one variable (cod_stomach_content_prop_large)
cod_prey_long_large %>% filter(cod_stomach_content_prop_large == "NaN")
#> filter: removed all rows (100%)
#> # A tibble: 0 × 3
#> # … with 3 variables: name <chr>, id <chr>,
#> # cod_stomach_content_prop_large <dbl>
cod_prey_long_small %>% filter(cod_stomach_content_prop_small == "NaN")
#> filter: removed 1,552 rows (99%), 16 rows remaining
#> # A tibble: 16 × 3
#> name id cod_stomach_content_prop_small
#> <chr> <chr> <dbl>
#> 1 amphipoda_tot 2017_1_43G9 NaN
#> 2 bivalvia_tot 2017_1_43G9 NaN
#> 3 clupea_harengus_tot 2017_1_43G9 NaN
#> 4 clupeidae_tot 2017_1_43G9 NaN
#> 5 gadiformes_tot 2017_1_43G9 NaN
#> 6 gadus_morhua_tot 2017_1_43G9 NaN
#> 7 gobiidae_tot 2017_1_43G9 NaN
#> 8 mysidae_tot 2017_1_43G9 NaN
#> 9 non_bio_tot 2017_1_43G9 NaN
#> 10 other_crustacea_tot 2017_1_43G9 NaN
#> 11 other_pisces_tot 2017_1_43G9 NaN
#> 12 other_tot 2017_1_43G9 NaN
#> 13 platichthys_flesus_tot 2017_1_43G9 NaN
#> 14 polychaeta_tot 2017_1_43G9 NaN
#> 15 saduria_entomon_tot 2017_1_43G9 NaN
#> 16 sprattus_sprattus_tot 2017_1_43G9 NaN
cod_prey_long_small <- cod_prey_long_small %>%
mutate(cod_stomach_content_prop_small = ifelse(cod_stomach_content_prop_small == "NaN",
0,
cod_stomach_content_prop_small))
#> mutate: changed 16 values (1%) of 'cod_stomach_content_prop_small' (16 fewer NA)
# Calculate Schoener index
schoener2 <- fle_prey_long %>%
left_join(cod_prey_long_small) %>%
left_join(cod_prey_long_large) %>%
drop_na(name) %>%
drop_na(fle_stomach_content_prop) %>%
drop_na(cod_stomach_content_prop_small) %>%
drop_na(cod_stomach_content_prop_large) %>%
group_by(year, quarter, ices_rect) %>%
summarise(schoener_f_sc = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop_small))),
schoener_f_lc = 1 - 0.5*(sum(abs(fle_stomach_content_prop - cod_stomach_content_prop_large))),
schoener_sc_lc = 1 - 0.5*(sum(abs(cod_stomach_content_prop_small - cod_stomach_content_prop_large)))) %>%
ungroup()
#> Joining, by = c("name", "id")
#> left_join: added one column (cod_stomach_content_prop_small)
#> > rows only in x 64
#> > rows only in y ( 336)
#> > matched rows 1,232
#> > =======
#> > rows total 1,296
#> Joining, by = c("name", "id")
#> left_join: added one column (cod_stomach_content_prop_large)
#> > rows only in x 48
#> > rows only in y ( 384)
#> > matched rows 1,248
#> > =======
#> > rows total 1,296
#> drop_na: no rows removed
#> drop_na: no rows removed
#> drop_na: removed 64 rows (5%), 1,232 rows remaining
#> drop_na: removed 16 rows (1%), 1,216 rows remaining
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 76 rows and 6 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
schoener_long2 <- schoener2 %>%
pivot_longer(4:6, names_to = "schoener_combination", values_to = "schoener")
#> pivot_longer: reorganized (schoener_f_sc, schoener_f_lc, schoener_sc_lc) into (schoener_combination, schoener) [was 76x6, now 228x5]
ggplot(schoener_long2, aes(schoener_combination, schoener, fill = factor(schoener_combination),
color = factor(schoener_combination))) +
ggdist::stat_halfeye(adjust = 0.5, justification = -0.1, .width = 0, point_colour = NA, alpha = 0.8,
show.legend = FALSE) +
geom_boxplot(width = 0.12, outlier.color = NA, alpha = 0.5, show.legend = FALSE) +
ggdist::stat_dots(side = "left", justification = 1.1, alpha = 0.6) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "",
labels = c("Flounder-Cod (L)", "Flounder-Cod (S)", "Cod (S)-Cod (L)")) +
guides(color = FALSE, fill = guide_legend(override.aes = list(
shape = 21, size = 2, fill = brewer.pal(n = 3, name = "Dark2"), color = brewer.pal(n = 3, name = "Dark2")))) +
coord_flip() +
labs(y = "Value", x = "") +
theme(legend.position = "bottom",
legend.direction = "horizontal",
axis.text.y = element_blank()) +
ggtitle("Schoeners's overlap index") +
NULL
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
ggsave("figures/schoener_size_groups_data.png", width = 6.5, height = 6.5, dpi = 600)
# Prepare for analysis
schoener_long2 <- schoener_long2 %>%
mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
schoener2 = ifelse(schoener == 1, 1-0.0001, schoener2),
year_f = as.factor(year),
quarter_f = as.factor(quarter),
ices_rect_f = as.factor(ices_rect),
schoener_combination_f = as.factor(schoener_combination))
#> mutate: new variable 'schoener2' (double) with 204 unique values and 0% NA
#> new variable 'year_f' (factor) with 6 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA
#> new variable 'schoener_combination_f' (factor) with 3 unique values and 0% NA
brms models of diversity and overlap indices# fit
m_schoen_beta <- brm(
bf(schoener2 ~ 0 + year_f + quarter_f + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f + predfle_density_tot_sc + predcod_density_tot_sc),
data = schoener, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_beta)
loo_m_schoen_beta <- loo(m_schoen_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
# Evaluate fit and convergence etc.
# PP check
pp_check(m_schoen_beta, ndraws = 50) +
theme_light(base_size = 20) +
scale_color_brewer(palette = "Dark2", name = "") +
NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/schoener_pp_check.png", width = 6.5, height = 6.5, dpi = 600)
# Chain convergence
posterior <- as.array(m_schoen_beta)
dimnames(posterior)
#> $iteration
#> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
#> [11] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
#> [21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
#> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
#> [41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50"
#> [51] "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
#> [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
#> [71] "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
#> [81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
#> [91] "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
#> [101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
#> [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
#> [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130"
#> [131] "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
#> [141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150"
#> [151] "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
#> [161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170"
#> [171] "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
#> [181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190"
#> [191] "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
#> [201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210"
#> [211] "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
#> [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230"
#> [231] "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
#> [241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250"
#> [251] "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
#> [261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270"
#> [271] "271" "272" "273" "274" "275" "276" "277" "278" "279" "280"
#> [281] "281" "282" "283" "284" "285" "286" "287" "288" "289" "290"
#> [291] "291" "292" "293" "294" "295" "296" "297" "298" "299" "300"
#> [301] "301" "302" "303" "304" "305" "306" "307" "308" "309" "310"
#> [311] "311" "312" "313" "314" "315" "316" "317" "318" "319" "320"
#> [321] "321" "322" "323" "324" "325" "326" "327" "328" "329" "330"
#> [331] "331" "332" "333" "334" "335" "336" "337" "338" "339" "340"
#> [341] "341" "342" "343" "344" "345" "346" "347" "348" "349" "350"
#> [351] "351" "352" "353" "354" "355" "356" "357" "358" "359" "360"
#> [361] "361" "362" "363" "364" "365" "366" "367" "368" "369" "370"
#> [371] "371" "372" "373" "374" "375" "376" "377" "378" "379" "380"
#> [381] "381" "382" "383" "384" "385" "386" "387" "388" "389" "390"
#> [391] "391" "392" "393" "394" "395" "396" "397" "398" "399" "400"
#> [401] "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
#> [411] "411" "412" "413" "414" "415" "416" "417" "418" "419" "420"
#> [421] "421" "422" "423" "424" "425" "426" "427" "428" "429" "430"
#> [431] "431" "432" "433" "434" "435" "436" "437" "438" "439" "440"
#> [441] "441" "442" "443" "444" "445" "446" "447" "448" "449" "450"
#> [451] "451" "452" "453" "454" "455" "456" "457" "458" "459" "460"
#> [461] "461" "462" "463" "464" "465" "466" "467" "468" "469" "470"
#> [471] "471" "472" "473" "474" "475" "476" "477" "478" "479" "480"
#> [481] "481" "482" "483" "484" "485" "486" "487" "488" "489" "490"
#> [491] "491" "492" "493" "494" "495" "496" "497" "498" "499" "500"
#> [501] "501" "502" "503" "504" "505" "506" "507" "508" "509" "510"
#> [511] "511" "512" "513" "514" "515" "516" "517" "518" "519" "520"
#> [521] "521" "522" "523" "524" "525" "526" "527" "528" "529" "530"
#> [531] "531" "532" "533" "534" "535" "536" "537" "538" "539" "540"
#> [541] "541" "542" "543" "544" "545" "546" "547" "548" "549" "550"
#> [551] "551" "552" "553" "554" "555" "556" "557" "558" "559" "560"
#> [561] "561" "562" "563" "564" "565" "566" "567" "568" "569" "570"
#> [571] "571" "572" "573" "574" "575" "576" "577" "578" "579" "580"
#> [581] "581" "582" "583" "584" "585" "586" "587" "588" "589" "590"
#> [591] "591" "592" "593" "594" "595" "596" "597" "598" "599" "600"
#> [601] "601" "602" "603" "604" "605" "606" "607" "608" "609" "610"
#> [611] "611" "612" "613" "614" "615" "616" "617" "618" "619" "620"
#> [621] "621" "622" "623" "624" "625" "626" "627" "628" "629" "630"
#> [631] "631" "632" "633" "634" "635" "636" "637" "638" "639" "640"
#> [641] "641" "642" "643" "644" "645" "646" "647" "648" "649" "650"
#> [651] "651" "652" "653" "654" "655" "656" "657" "658" "659" "660"
#> [661] "661" "662" "663" "664" "665" "666" "667" "668" "669" "670"
#> [671] "671" "672" "673" "674" "675" "676" "677" "678" "679" "680"
#> [681] "681" "682" "683" "684" "685" "686" "687" "688" "689" "690"
#> [691] "691" "692" "693" "694" "695" "696" "697" "698" "699" "700"
#> [701] "701" "702" "703" "704" "705" "706" "707" "708" "709" "710"
#> [711] "711" "712" "713" "714" "715" "716" "717" "718" "719" "720"
#> [721] "721" "722" "723" "724" "725" "726" "727" "728" "729" "730"
#> [731] "731" "732" "733" "734" "735" "736" "737" "738" "739" "740"
#> [741] "741" "742" "743" "744" "745" "746" "747" "748" "749" "750"
#> [751] "751" "752" "753" "754" "755" "756" "757" "758" "759" "760"
#> [761] "761" "762" "763" "764" "765" "766" "767" "768" "769" "770"
#> [771] "771" "772" "773" "774" "775" "776" "777" "778" "779" "780"
#> [781] "781" "782" "783" "784" "785" "786" "787" "788" "789" "790"
#> [791] "791" "792" "793" "794" "795" "796" "797" "798" "799" "800"
#> [801] "801" "802" "803" "804" "805" "806" "807" "808" "809" "810"
#> [811] "811" "812" "813" "814" "815" "816" "817" "818" "819" "820"
#> [821] "821" "822" "823" "824" "825" "826" "827" "828" "829" "830"
#> [831] "831" "832" "833" "834" "835" "836" "837" "838" "839" "840"
#> [841] "841" "842" "843" "844" "845" "846" "847" "848" "849" "850"
#> [851] "851" "852" "853" "854" "855" "856" "857" "858" "859" "860"
#> [861] "861" "862" "863" "864" "865" "866" "867" "868" "869" "870"
#> [871] "871" "872" "873" "874" "875" "876" "877" "878" "879" "880"
#> [881] "881" "882" "883" "884" "885" "886" "887" "888" "889" "890"
#> [891] "891" "892" "893" "894" "895" "896" "897" "898" "899" "900"
#> [901] "901" "902" "903" "904" "905" "906" "907" "908" "909" "910"
#> [911] "911" "912" "913" "914" "915" "916" "917" "918" "919" "920"
#> [921] "921" "922" "923" "924" "925" "926" "927" "928" "929" "930"
#> [931] "931" "932" "933" "934" "935" "936" "937" "938" "939" "940"
#> [941] "941" "942" "943" "944" "945" "946" "947" "948" "949" "950"
#> [951] "951" "952" "953" "954" "955" "956" "957" "958" "959" "960"
#> [961] "961" "962" "963" "964" "965" "966" "967" "968" "969" "970"
#> [971] "971" "972" "973" "974" "975" "976" "977" "978" "979" "980"
#> [981] "981" "982" "983" "984" "985" "986" "987" "988" "989" "990"
#> [991] "991" "992" "993" "994" "995" "996" "997" "998" "999" "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#>
#> $chain
#> [1] "1" "2" "3" "4"
#>
#> $variable
#> [1] "b_year_f2015" "b_year_f2016"
#> [3] "b_year_f2017" "b_year_f2018"
#> [5] "b_year_f2019" "b_year_f2020"
#> [7] "b_quarter_f4" "b_predfle_density_tot_sc"
#> [9] "b_predcod_density_tot_sc" "b_phi_year_f2015"
#> [11] "b_phi_year_f2016" "b_phi_year_f2017"
#> [13] "b_phi_year_f2018" "b_phi_year_f2019"
#> [15] "b_phi_year_f2020" "b_phi_quarter_f4"
#> [17] "b_phi_predfle_density_tot_sc" "b_phi_predcod_density_tot_sc"
#> [19] "sd_ices_rect_f__Intercept" "r_ices_rect_f[39G4,Intercept]"
#> [21] "r_ices_rect_f[40G4,Intercept]" "r_ices_rect_f[40G5,Intercept]"
#> [23] "r_ices_rect_f[40G6,Intercept]" "r_ices_rect_f[40G7,Intercept]"
#> [25] "r_ices_rect_f[40G8,Intercept]" "r_ices_rect_f[41G4,Intercept]"
#> [27] "r_ices_rect_f[41G6,Intercept]" "r_ices_rect_f[41G7,Intercept]"
#> [29] "r_ices_rect_f[41G8,Intercept]" "r_ices_rect_f[41G9,Intercept]"
#> [31] "r_ices_rect_f[42G6,Intercept]" "r_ices_rect_f[42G7,Intercept]"
#> [33] "r_ices_rect_f[43G6,Intercept]" "r_ices_rect_f[43G7,Intercept]"
#> [35] "r_ices_rect_f[43G8,Intercept]" "r_ices_rect_f[43G9,Intercept]"
#> [37] "r_ices_rect_f[44G9,Intercept]" "lp__"
#> [39] "z_1[1,1]" "z_1[1,2]"
#> [41] "z_1[1,3]" "z_1[1,4]"
#> [43] "z_1[1,5]" "z_1[1,6]"
#> [45] "z_1[1,7]" "z_1[1,8]"
#> [47] "z_1[1,9]" "z_1[1,10]"
#> [49] "z_1[1,11]" "z_1[1,12]"
#> [51] "z_1[1,13]" "z_1[1,14]"
#> [53] "z_1[1,15]" "z_1[1,16]"
#> [55] "z_1[1,17]" "z_1[1,18]"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))
mcmc_trace(posterior,
pars = c("b_year_f2015", "b_year_f2016", "b_year_f2017", "b_year_f2018",
"b_year_f2018", "b_year_f2019", "b_year_f2020",
"b_predfle_density_tot_sc", "b_predcod_density_tot_sc",
"b_quarter_f4", "b_phi_quarter_f4",
"b_phi_year_f2015", "b_phi_year_f2016", "b_phi_year_f2017",
"b_phi_year_f2018", "b_phi_year_f2019", "b_phi_year_f2020",
"b_phi_predfle_density_tot_sc", "b_phi_predcod_density_tot_sc",
"sd_ices_rect_f__Intercept"),
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12), strip.text = element_text(size = 4),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/schoener_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Working with the posterior
posterior_m_schoen_beta <- m_schoen_beta %>%
gather_draws(`b_.*`, regex = TRUE)
ggplot(posterior_m_schoen_beta, aes(x = .value, y = fct_rev(.variable))) +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.7) +
stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95), point_interval = "median_hdi") +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Variable") +
NULL
# Marginal effects of densities
m_schoen_beta_pred_cod <- m_schoen_beta %>%
epred_draws(newdata = tibble(predcod_density_tot_sc =
seq(min(schoener$predcod_density_tot_sc),
max(schoener$predcod_density_tot_sc),
length.out = 1000),
predfle_density_tot_sc = 0,
year_f = factor(2018),
quarter_f = factor(1)),
re_formula = NA)
m_schoen_beta_pred_fle <- m_schoen_beta %>%
epred_draws(newdata = tibble(predfle_density_tot_sc =
seq(min(schoener$predfle_density_tot_sc),
max(schoener$predfle_density_tot_sc),
length.out = 1000),
predcod_density_tot_sc = 0,
year_f = factor(2018),
quarter_f = factor(1)),
re_formula = NA)
p1 <- ggplot(m_schoen_beta_pred_fle, aes(x = predfle_density_tot_sc, y = .epred)) +
stat_lineribbon() +
scale_fill_brewer(palette = "Purples") +
theme_light(base_size = 14) +
theme(legend.position = "bottom") +
guides(fill = FALSE) +
labs(x = "Scaled flounder density", y = NULL
#, caption = "80% and 95% credible intervals shown in black"
) +
NULL
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
p2 <- ggplot(m_schoen_beta_pred_cod, aes(x = predcod_density_tot_sc, y = .epred)) +
stat_lineribbon() +
scale_fill_brewer(palette = "Purples", name = "Credible interval") +
theme_light(base_size = 14) +
theme(legend.position = "bottom") +
guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
labs(x = "Scaled cod density", y = NULL
, caption = "80% and 95% credible intervals shown in black"
) +
NULL
p1/p2
ggsave("figures/schoener_dens_marginal.png", width = 6.5, height = 6.5, dpi = 600)
m_schoen_size_beta <- brm(
bf(schoener2 ~ 0 + schoener_combination_f + year_f + quarter_f + (1|ices_rect_f),
phi ~ 0 + schoener_combination_f + year_f + quarter_f),
data = schoener_long2, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_size_beta)
loo_m_schoen_size_beta <- loo(m_schoen_size_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
# Evaluate fit and convergence etc.
# PP check
pp_check(m_schoen_size_beta, ndraws = 50) +
theme_light(base_size = 20) +
scale_color_brewer(palette = "Dark2", name = "") +
NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/schoener_size_pp_check.png", width = 6.5, height = 6.5, dpi = 600)
# Chain convergence
posterior <- as.array(m_schoen_size_beta)
dimnames(posterior)
#> $iteration
#> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
#> [11] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
#> [21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
#> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
#> [41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50"
#> [51] "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
#> [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
#> [71] "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
#> [81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
#> [91] "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
#> [101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
#> [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
#> [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130"
#> [131] "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
#> [141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150"
#> [151] "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
#> [161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170"
#> [171] "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
#> [181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190"
#> [191] "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
#> [201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210"
#> [211] "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
#> [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230"
#> [231] "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
#> [241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250"
#> [251] "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
#> [261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270"
#> [271] "271" "272" "273" "274" "275" "276" "277" "278" "279" "280"
#> [281] "281" "282" "283" "284" "285" "286" "287" "288" "289" "290"
#> [291] "291" "292" "293" "294" "295" "296" "297" "298" "299" "300"
#> [301] "301" "302" "303" "304" "305" "306" "307" "308" "309" "310"
#> [311] "311" "312" "313" "314" "315" "316" "317" "318" "319" "320"
#> [321] "321" "322" "323" "324" "325" "326" "327" "328" "329" "330"
#> [331] "331" "332" "333" "334" "335" "336" "337" "338" "339" "340"
#> [341] "341" "342" "343" "344" "345" "346" "347" "348" "349" "350"
#> [351] "351" "352" "353" "354" "355" "356" "357" "358" "359" "360"
#> [361] "361" "362" "363" "364" "365" "366" "367" "368" "369" "370"
#> [371] "371" "372" "373" "374" "375" "376" "377" "378" "379" "380"
#> [381] "381" "382" "383" "384" "385" "386" "387" "388" "389" "390"
#> [391] "391" "392" "393" "394" "395" "396" "397" "398" "399" "400"
#> [401] "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
#> [411] "411" "412" "413" "414" "415" "416" "417" "418" "419" "420"
#> [421] "421" "422" "423" "424" "425" "426" "427" "428" "429" "430"
#> [431] "431" "432" "433" "434" "435" "436" "437" "438" "439" "440"
#> [441] "441" "442" "443" "444" "445" "446" "447" "448" "449" "450"
#> [451] "451" "452" "453" "454" "455" "456" "457" "458" "459" "460"
#> [461] "461" "462" "463" "464" "465" "466" "467" "468" "469" "470"
#> [471] "471" "472" "473" "474" "475" "476" "477" "478" "479" "480"
#> [481] "481" "482" "483" "484" "485" "486" "487" "488" "489" "490"
#> [491] "491" "492" "493" "494" "495" "496" "497" "498" "499" "500"
#> [501] "501" "502" "503" "504" "505" "506" "507" "508" "509" "510"
#> [511] "511" "512" "513" "514" "515" "516" "517" "518" "519" "520"
#> [521] "521" "522" "523" "524" "525" "526" "527" "528" "529" "530"
#> [531] "531" "532" "533" "534" "535" "536" "537" "538" "539" "540"
#> [541] "541" "542" "543" "544" "545" "546" "547" "548" "549" "550"
#> [551] "551" "552" "553" "554" "555" "556" "557" "558" "559" "560"
#> [561] "561" "562" "563" "564" "565" "566" "567" "568" "569" "570"
#> [571] "571" "572" "573" "574" "575" "576" "577" "578" "579" "580"
#> [581] "581" "582" "583" "584" "585" "586" "587" "588" "589" "590"
#> [591] "591" "592" "593" "594" "595" "596" "597" "598" "599" "600"
#> [601] "601" "602" "603" "604" "605" "606" "607" "608" "609" "610"
#> [611] "611" "612" "613" "614" "615" "616" "617" "618" "619" "620"
#> [621] "621" "622" "623" "624" "625" "626" "627" "628" "629" "630"
#> [631] "631" "632" "633" "634" "635" "636" "637" "638" "639" "640"
#> [641] "641" "642" "643" "644" "645" "646" "647" "648" "649" "650"
#> [651] "651" "652" "653" "654" "655" "656" "657" "658" "659" "660"
#> [661] "661" "662" "663" "664" "665" "666" "667" "668" "669" "670"
#> [671] "671" "672" "673" "674" "675" "676" "677" "678" "679" "680"
#> [681] "681" "682" "683" "684" "685" "686" "687" "688" "689" "690"
#> [691] "691" "692" "693" "694" "695" "696" "697" "698" "699" "700"
#> [701] "701" "702" "703" "704" "705" "706" "707" "708" "709" "710"
#> [711] "711" "712" "713" "714" "715" "716" "717" "718" "719" "720"
#> [721] "721" "722" "723" "724" "725" "726" "727" "728" "729" "730"
#> [731] "731" "732" "733" "734" "735" "736" "737" "738" "739" "740"
#> [741] "741" "742" "743" "744" "745" "746" "747" "748" "749" "750"
#> [751] "751" "752" "753" "754" "755" "756" "757" "758" "759" "760"
#> [761] "761" "762" "763" "764" "765" "766" "767" "768" "769" "770"
#> [771] "771" "772" "773" "774" "775" "776" "777" "778" "779" "780"
#> [781] "781" "782" "783" "784" "785" "786" "787" "788" "789" "790"
#> [791] "791" "792" "793" "794" "795" "796" "797" "798" "799" "800"
#> [801] "801" "802" "803" "804" "805" "806" "807" "808" "809" "810"
#> [811] "811" "812" "813" "814" "815" "816" "817" "818" "819" "820"
#> [821] "821" "822" "823" "824" "825" "826" "827" "828" "829" "830"
#> [831] "831" "832" "833" "834" "835" "836" "837" "838" "839" "840"
#> [841] "841" "842" "843" "844" "845" "846" "847" "848" "849" "850"
#> [851] "851" "852" "853" "854" "855" "856" "857" "858" "859" "860"
#> [861] "861" "862" "863" "864" "865" "866" "867" "868" "869" "870"
#> [871] "871" "872" "873" "874" "875" "876" "877" "878" "879" "880"
#> [881] "881" "882" "883" "884" "885" "886" "887" "888" "889" "890"
#> [891] "891" "892" "893" "894" "895" "896" "897" "898" "899" "900"
#> [901] "901" "902" "903" "904" "905" "906" "907" "908" "909" "910"
#> [911] "911" "912" "913" "914" "915" "916" "917" "918" "919" "920"
#> [921] "921" "922" "923" "924" "925" "926" "927" "928" "929" "930"
#> [931] "931" "932" "933" "934" "935" "936" "937" "938" "939" "940"
#> [941] "941" "942" "943" "944" "945" "946" "947" "948" "949" "950"
#> [951] "951" "952" "953" "954" "955" "956" "957" "958" "959" "960"
#> [961] "961" "962" "963" "964" "965" "966" "967" "968" "969" "970"
#> [971] "971" "972" "973" "974" "975" "976" "977" "978" "979" "980"
#> [981] "981" "982" "983" "984" "985" "986" "987" "988" "989" "990"
#> [991] "991" "992" "993" "994" "995" "996" "997" "998" "999" "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#>
#> $chain
#> [1] "1" "2" "3" "4"
#>
#> $variable
#> [1] "b_schoener_combination_fschoener_f_lc"
#> [2] "b_schoener_combination_fschoener_f_sc"
#> [3] "b_schoener_combination_fschoener_sc_lc"
#> [4] "b_year_f2016"
#> [5] "b_year_f2017"
#> [6] "b_year_f2018"
#> [7] "b_year_f2019"
#> [8] "b_year_f2020"
#> [9] "b_quarter_f4"
#> [10] "b_phi_schoener_combination_fschoener_f_lc"
#> [11] "b_phi_schoener_combination_fschoener_f_sc"
#> [12] "b_phi_schoener_combination_fschoener_sc_lc"
#> [13] "b_phi_year_f2016"
#> [14] "b_phi_year_f2017"
#> [15] "b_phi_year_f2018"
#> [16] "b_phi_year_f2019"
#> [17] "b_phi_year_f2020"
#> [18] "b_phi_quarter_f4"
#> [19] "sd_ices_rect_f__Intercept"
#> [20] "r_ices_rect_f[39G4,Intercept]"
#> [21] "r_ices_rect_f[40G4,Intercept]"
#> [22] "r_ices_rect_f[40G5,Intercept]"
#> [23] "r_ices_rect_f[40G6,Intercept]"
#> [24] "r_ices_rect_f[40G7,Intercept]"
#> [25] "r_ices_rect_f[40G8,Intercept]"
#> [26] "r_ices_rect_f[41G4,Intercept]"
#> [27] "r_ices_rect_f[41G6,Intercept]"
#> [28] "r_ices_rect_f[41G7,Intercept]"
#> [29] "r_ices_rect_f[41G8,Intercept]"
#> [30] "r_ices_rect_f[41G9,Intercept]"
#> [31] "r_ices_rect_f[42G6,Intercept]"
#> [32] "r_ices_rect_f[42G7,Intercept]"
#> [33] "r_ices_rect_f[43G6,Intercept]"
#> [34] "r_ices_rect_f[43G7,Intercept]"
#> [35] "r_ices_rect_f[43G8,Intercept]"
#> [36] "r_ices_rect_f[43G9,Intercept]"
#> [37] "r_ices_rect_f[44G9,Intercept]"
#> [38] "lp__"
#> [39] "z_1[1,1]"
#> [40] "z_1[1,2]"
#> [41] "z_1[1,3]"
#> [42] "z_1[1,4]"
#> [43] "z_1[1,5]"
#> [44] "z_1[1,6]"
#> [45] "z_1[1,7]"
#> [46] "z_1[1,8]"
#> [47] "z_1[1,9]"
#> [48] "z_1[1,10]"
#> [49] "z_1[1,11]"
#> [50] "z_1[1,12]"
#> [51] "z_1[1,13]"
#> [52] "z_1[1,14]"
#> [53] "z_1[1,15]"
#> [54] "z_1[1,16]"
#> [55] "z_1[1,17]"
#> [56] "z_1[1,18]"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))
mcmc_trace(posterior,
pars = c("b_schoener_combination_fschoener_f_lc",
"b_schoener_combination_fschoener_f_sc",
"b_schoener_combination_fschoener_sc_lc",
"b_year_f2016", "b_year_f2017", "b_year_f2018", "b_year_f2019",
"b_year_f2020", "b_quarter_f4", "b_phi_quarter_f4",
"b_phi_year_f2016", "b_phi_year_f2017", "b_phi_year_f2018",
"b_phi_year_f2019", "b_phi_year_f2020"),
facet_args = list(ncol = 2, strip.position = "left")) +
theme(text = element_text(size = 12), strip.text = element_text(size = 4),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/schoener_size_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Working with the posterior
posterior_m_schoen_size_beta <- m_schoen_size_beta %>%
gather_draws(`b_.*`, regex = TRUE)
ggplot(posterior_m_schoen_size_beta, aes(x = .value, y = fct_rev(.variable))) +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.7) +
stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95), point_interval = "median_hdi") +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Variable") +
NULL
# Marginal effects of Schoener combination variable
m_schoen_size_beta_pred <- m_schoen_size_beta %>%
epred_draws(newdata = tibble(schoener_combination_f = c("schoener_f_sc", "schoener_f_lc", "schoener_sc_lc"),
year_f = factor(2018),
quarter_f = factor(1)),
re_formula = NA) %>%
mutate(schoener_combination_f = ifelse(schoener_combination_f == "schoener_f_sc",
"Flounder-Cod (S)", schoener_combination_f),
schoener_combination_f = ifelse(schoener_combination_f == "schoener_f_lc",
"Flounder-Cod (L)", schoener_combination_f),
schoener_combination_f = ifelse(schoener_combination_f == "schoener_sc_lc",
"Cod (S)-Cod (L)", schoener_combination_f))
#> mutate (grouped): changed 24,000 values (100%) of 'schoener_combination_f' (0 new NA)
ggplot(m_schoen_size_beta_pred, aes(x = .epred, fill = schoener_combination_f)) +
stat_halfeye(.width = c(0.75, 0.95), point_interval = "median_hdi", alpha = 0.5,
position = position_dodge(width = 0.01)) +
coord_cartesian(ylim = c(0, 0.3)) +
scale_fill_brewer(palette = "Dark2", name = "Species-size combination") +
theme_light(base_size = 14) +
theme(legend.position = "top") +
geom_vline(xintercept = 0.6, linetype = 2, color = "gray20") +
guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
labs(x = "Predicted Schoener's overlap index", y = NULL
, caption = "80% and 95% credible intervals shown in black"
) +
NULL
ggsave("figures/schoener_size_marginal.png", width = 6.5, height = 6.5, dpi = 600)
# Cod model
m_levin_beta_cod <- brm(
bf(levin_cod2 ~ 1,
phi ~ 1),
data = levin, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta_cod)
loo_m_levin_beta_cod <- loo(m_levin_beta_cod)
# Flounder models
m_levin_beta_fle <- brm(
bf(levin_fle2 ~ 1,
phi ~ 1),
data = levin, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> recompiling to avoid crashing R session
#> Start sampling
plot(m_levin_beta_fle)
loo_m_levin_beta_fle <- loo(m_levin_beta_fle)
# Pooled model of Levin index and species as covariate
colnames(levin)
#> [1] "year" "quarter" "area" "levin_cod" "levin_fle"
#> [6] "levin_cod2" "levin_fle2" "year_f" "quarter_f" "area_f"
levin2 <- levin %>% pivot_longer(4:5, names_to = "species", values_to = "levins")
#> pivot_longer: reorganized (levin_cod, levin_fle) into (species, levins) [was 22x10, now 44x10]
colnames(levin2)
#> [1] "year" "quarter" "area" "levin_cod2" "levin_fle2"
#> [6] "year_f" "quarter_f" "area_f" "species" "levins"
m_levin_beta <- brm(
bf(levins ~ 0 + species,
phi ~ 0 + species),
data = levin2, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
chains = 4, iter = 4000, cores = 4, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta)
loo_m_levin_beta <- loo(m_levin_beta)
# Evaluate fit and convergence etc.
# PP check
pp_check(m_levin_beta, ndraws = 50) +
theme_light(base_size = 20) +
scale_color_brewer(palette = "Dark2", name = "") +
NULL
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/levin_pp_check.png", width = 6.5, height = 6.5, dpi = 600)
# Chain convergence
posterior <- as.array(m_levin_beta)
dimnames(posterior)
#> $iteration
#> [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
#> [11] "11" "12" "13" "14" "15" "16" "17" "18" "19" "20"
#> [21] "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
#> [31] "31" "32" "33" "34" "35" "36" "37" "38" "39" "40"
#> [41] "41" "42" "43" "44" "45" "46" "47" "48" "49" "50"
#> [51] "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
#> [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70"
#> [71] "71" "72" "73" "74" "75" "76" "77" "78" "79" "80"
#> [81] "81" "82" "83" "84" "85" "86" "87" "88" "89" "90"
#> [91] "91" "92" "93" "94" "95" "96" "97" "98" "99" "100"
#> [101] "101" "102" "103" "104" "105" "106" "107" "108" "109" "110"
#> [111] "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
#> [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130"
#> [131] "131" "132" "133" "134" "135" "136" "137" "138" "139" "140"
#> [141] "141" "142" "143" "144" "145" "146" "147" "148" "149" "150"
#> [151] "151" "152" "153" "154" "155" "156" "157" "158" "159" "160"
#> [161] "161" "162" "163" "164" "165" "166" "167" "168" "169" "170"
#> [171] "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
#> [181] "181" "182" "183" "184" "185" "186" "187" "188" "189" "190"
#> [191] "191" "192" "193" "194" "195" "196" "197" "198" "199" "200"
#> [201] "201" "202" "203" "204" "205" "206" "207" "208" "209" "210"
#> [211] "211" "212" "213" "214" "215" "216" "217" "218" "219" "220"
#> [221] "221" "222" "223" "224" "225" "226" "227" "228" "229" "230"
#> [231] "231" "232" "233" "234" "235" "236" "237" "238" "239" "240"
#> [241] "241" "242" "243" "244" "245" "246" "247" "248" "249" "250"
#> [251] "251" "252" "253" "254" "255" "256" "257" "258" "259" "260"
#> [261] "261" "262" "263" "264" "265" "266" "267" "268" "269" "270"
#> [271] "271" "272" "273" "274" "275" "276" "277" "278" "279" "280"
#> [281] "281" "282" "283" "284" "285" "286" "287" "288" "289" "290"
#> [291] "291" "292" "293" "294" "295" "296" "297" "298" "299" "300"
#> [301] "301" "302" "303" "304" "305" "306" "307" "308" "309" "310"
#> [311] "311" "312" "313" "314" "315" "316" "317" "318" "319" "320"
#> [321] "321" "322" "323" "324" "325" "326" "327" "328" "329" "330"
#> [331] "331" "332" "333" "334" "335" "336" "337" "338" "339" "340"
#> [341] "341" "342" "343" "344" "345" "346" "347" "348" "349" "350"
#> [351] "351" "352" "353" "354" "355" "356" "357" "358" "359" "360"
#> [361] "361" "362" "363" "364" "365" "366" "367" "368" "369" "370"
#> [371] "371" "372" "373" "374" "375" "376" "377" "378" "379" "380"
#> [381] "381" "382" "383" "384" "385" "386" "387" "388" "389" "390"
#> [391] "391" "392" "393" "394" "395" "396" "397" "398" "399" "400"
#> [401] "401" "402" "403" "404" "405" "406" "407" "408" "409" "410"
#> [411] "411" "412" "413" "414" "415" "416" "417" "418" "419" "420"
#> [421] "421" "422" "423" "424" "425" "426" "427" "428" "429" "430"
#> [431] "431" "432" "433" "434" "435" "436" "437" "438" "439" "440"
#> [441] "441" "442" "443" "444" "445" "446" "447" "448" "449" "450"
#> [451] "451" "452" "453" "454" "455" "456" "457" "458" "459" "460"
#> [461] "461" "462" "463" "464" "465" "466" "467" "468" "469" "470"
#> [471] "471" "472" "473" "474" "475" "476" "477" "478" "479" "480"
#> [481] "481" "482" "483" "484" "485" "486" "487" "488" "489" "490"
#> [491] "491" "492" "493" "494" "495" "496" "497" "498" "499" "500"
#> [501] "501" "502" "503" "504" "505" "506" "507" "508" "509" "510"
#> [511] "511" "512" "513" "514" "515" "516" "517" "518" "519" "520"
#> [521] "521" "522" "523" "524" "525" "526" "527" "528" "529" "530"
#> [531] "531" "532" "533" "534" "535" "536" "537" "538" "539" "540"
#> [541] "541" "542" "543" "544" "545" "546" "547" "548" "549" "550"
#> [551] "551" "552" "553" "554" "555" "556" "557" "558" "559" "560"
#> [561] "561" "562" "563" "564" "565" "566" "567" "568" "569" "570"
#> [571] "571" "572" "573" "574" "575" "576" "577" "578" "579" "580"
#> [581] "581" "582" "583" "584" "585" "586" "587" "588" "589" "590"
#> [591] "591" "592" "593" "594" "595" "596" "597" "598" "599" "600"
#> [601] "601" "602" "603" "604" "605" "606" "607" "608" "609" "610"
#> [611] "611" "612" "613" "614" "615" "616" "617" "618" "619" "620"
#> [621] "621" "622" "623" "624" "625" "626" "627" "628" "629" "630"
#> [631] "631" "632" "633" "634" "635" "636" "637" "638" "639" "640"
#> [641] "641" "642" "643" "644" "645" "646" "647" "648" "649" "650"
#> [651] "651" "652" "653" "654" "655" "656" "657" "658" "659" "660"
#> [661] "661" "662" "663" "664" "665" "666" "667" "668" "669" "670"
#> [671] "671" "672" "673" "674" "675" "676" "677" "678" "679" "680"
#> [681] "681" "682" "683" "684" "685" "686" "687" "688" "689" "690"
#> [691] "691" "692" "693" "694" "695" "696" "697" "698" "699" "700"
#> [701] "701" "702" "703" "704" "705" "706" "707" "708" "709" "710"
#> [711] "711" "712" "713" "714" "715" "716" "717" "718" "719" "720"
#> [721] "721" "722" "723" "724" "725" "726" "727" "728" "729" "730"
#> [731] "731" "732" "733" "734" "735" "736" "737" "738" "739" "740"
#> [741] "741" "742" "743" "744" "745" "746" "747" "748" "749" "750"
#> [751] "751" "752" "753" "754" "755" "756" "757" "758" "759" "760"
#> [761] "761" "762" "763" "764" "765" "766" "767" "768" "769" "770"
#> [771] "771" "772" "773" "774" "775" "776" "777" "778" "779" "780"
#> [781] "781" "782" "783" "784" "785" "786" "787" "788" "789" "790"
#> [791] "791" "792" "793" "794" "795" "796" "797" "798" "799" "800"
#> [801] "801" "802" "803" "804" "805" "806" "807" "808" "809" "810"
#> [811] "811" "812" "813" "814" "815" "816" "817" "818" "819" "820"
#> [821] "821" "822" "823" "824" "825" "826" "827" "828" "829" "830"
#> [831] "831" "832" "833" "834" "835" "836" "837" "838" "839" "840"
#> [841] "841" "842" "843" "844" "845" "846" "847" "848" "849" "850"
#> [851] "851" "852" "853" "854" "855" "856" "857" "858" "859" "860"
#> [861] "861" "862" "863" "864" "865" "866" "867" "868" "869" "870"
#> [871] "871" "872" "873" "874" "875" "876" "877" "878" "879" "880"
#> [881] "881" "882" "883" "884" "885" "886" "887" "888" "889" "890"
#> [891] "891" "892" "893" "894" "895" "896" "897" "898" "899" "900"
#> [901] "901" "902" "903" "904" "905" "906" "907" "908" "909" "910"
#> [911] "911" "912" "913" "914" "915" "916" "917" "918" "919" "920"
#> [921] "921" "922" "923" "924" "925" "926" "927" "928" "929" "930"
#> [931] "931" "932" "933" "934" "935" "936" "937" "938" "939" "940"
#> [941] "941" "942" "943" "944" "945" "946" "947" "948" "949" "950"
#> [951] "951" "952" "953" "954" "955" "956" "957" "958" "959" "960"
#> [961] "961" "962" "963" "964" "965" "966" "967" "968" "969" "970"
#> [971] "971" "972" "973" "974" "975" "976" "977" "978" "979" "980"
#> [981] "981" "982" "983" "984" "985" "986" "987" "988" "989" "990"
#> [991] "991" "992" "993" "994" "995" "996" "997" "998" "999" "1000"
#> [1001] "1001" "1002" "1003" "1004" "1005" "1006" "1007" "1008" "1009" "1010"
#> [1011] "1011" "1012" "1013" "1014" "1015" "1016" "1017" "1018" "1019" "1020"
#> [1021] "1021" "1022" "1023" "1024" "1025" "1026" "1027" "1028" "1029" "1030"
#> [1031] "1031" "1032" "1033" "1034" "1035" "1036" "1037" "1038" "1039" "1040"
#> [1041] "1041" "1042" "1043" "1044" "1045" "1046" "1047" "1048" "1049" "1050"
#> [1051] "1051" "1052" "1053" "1054" "1055" "1056" "1057" "1058" "1059" "1060"
#> [1061] "1061" "1062" "1063" "1064" "1065" "1066" "1067" "1068" "1069" "1070"
#> [1071] "1071" "1072" "1073" "1074" "1075" "1076" "1077" "1078" "1079" "1080"
#> [1081] "1081" "1082" "1083" "1084" "1085" "1086" "1087" "1088" "1089" "1090"
#> [1091] "1091" "1092" "1093" "1094" "1095" "1096" "1097" "1098" "1099" "1100"
#> [1101] "1101" "1102" "1103" "1104" "1105" "1106" "1107" "1108" "1109" "1110"
#> [1111] "1111" "1112" "1113" "1114" "1115" "1116" "1117" "1118" "1119" "1120"
#> [1121] "1121" "1122" "1123" "1124" "1125" "1126" "1127" "1128" "1129" "1130"
#> [1131] "1131" "1132" "1133" "1134" "1135" "1136" "1137" "1138" "1139" "1140"
#> [1141] "1141" "1142" "1143" "1144" "1145" "1146" "1147" "1148" "1149" "1150"
#> [1151] "1151" "1152" "1153" "1154" "1155" "1156" "1157" "1158" "1159" "1160"
#> [1161] "1161" "1162" "1163" "1164" "1165" "1166" "1167" "1168" "1169" "1170"
#> [1171] "1171" "1172" "1173" "1174" "1175" "1176" "1177" "1178" "1179" "1180"
#> [1181] "1181" "1182" "1183" "1184" "1185" "1186" "1187" "1188" "1189" "1190"
#> [1191] "1191" "1192" "1193" "1194" "1195" "1196" "1197" "1198" "1199" "1200"
#> [1201] "1201" "1202" "1203" "1204" "1205" "1206" "1207" "1208" "1209" "1210"
#> [1211] "1211" "1212" "1213" "1214" "1215" "1216" "1217" "1218" "1219" "1220"
#> [1221] "1221" "1222" "1223" "1224" "1225" "1226" "1227" "1228" "1229" "1230"
#> [1231] "1231" "1232" "1233" "1234" "1235" "1236" "1237" "1238" "1239" "1240"
#> [1241] "1241" "1242" "1243" "1244" "1245" "1246" "1247" "1248" "1249" "1250"
#> [1251] "1251" "1252" "1253" "1254" "1255" "1256" "1257" "1258" "1259" "1260"
#> [1261] "1261" "1262" "1263" "1264" "1265" "1266" "1267" "1268" "1269" "1270"
#> [1271] "1271" "1272" "1273" "1274" "1275" "1276" "1277" "1278" "1279" "1280"
#> [1281] "1281" "1282" "1283" "1284" "1285" "1286" "1287" "1288" "1289" "1290"
#> [1291] "1291" "1292" "1293" "1294" "1295" "1296" "1297" "1298" "1299" "1300"
#> [1301] "1301" "1302" "1303" "1304" "1305" "1306" "1307" "1308" "1309" "1310"
#> [1311] "1311" "1312" "1313" "1314" "1315" "1316" "1317" "1318" "1319" "1320"
#> [1321] "1321" "1322" "1323" "1324" "1325" "1326" "1327" "1328" "1329" "1330"
#> [1331] "1331" "1332" "1333" "1334" "1335" "1336" "1337" "1338" "1339" "1340"
#> [1341] "1341" "1342" "1343" "1344" "1345" "1346" "1347" "1348" "1349" "1350"
#> [1351] "1351" "1352" "1353" "1354" "1355" "1356" "1357" "1358" "1359" "1360"
#> [1361] "1361" "1362" "1363" "1364" "1365" "1366" "1367" "1368" "1369" "1370"
#> [1371] "1371" "1372" "1373" "1374" "1375" "1376" "1377" "1378" "1379" "1380"
#> [1381] "1381" "1382" "1383" "1384" "1385" "1386" "1387" "1388" "1389" "1390"
#> [1391] "1391" "1392" "1393" "1394" "1395" "1396" "1397" "1398" "1399" "1400"
#> [1401] "1401" "1402" "1403" "1404" "1405" "1406" "1407" "1408" "1409" "1410"
#> [1411] "1411" "1412" "1413" "1414" "1415" "1416" "1417" "1418" "1419" "1420"
#> [1421] "1421" "1422" "1423" "1424" "1425" "1426" "1427" "1428" "1429" "1430"
#> [1431] "1431" "1432" "1433" "1434" "1435" "1436" "1437" "1438" "1439" "1440"
#> [1441] "1441" "1442" "1443" "1444" "1445" "1446" "1447" "1448" "1449" "1450"
#> [1451] "1451" "1452" "1453" "1454" "1455" "1456" "1457" "1458" "1459" "1460"
#> [1461] "1461" "1462" "1463" "1464" "1465" "1466" "1467" "1468" "1469" "1470"
#> [1471] "1471" "1472" "1473" "1474" "1475" "1476" "1477" "1478" "1479" "1480"
#> [1481] "1481" "1482" "1483" "1484" "1485" "1486" "1487" "1488" "1489" "1490"
#> [1491] "1491" "1492" "1493" "1494" "1495" "1496" "1497" "1498" "1499" "1500"
#> [1501] "1501" "1502" "1503" "1504" "1505" "1506" "1507" "1508" "1509" "1510"
#> [1511] "1511" "1512" "1513" "1514" "1515" "1516" "1517" "1518" "1519" "1520"
#> [1521] "1521" "1522" "1523" "1524" "1525" "1526" "1527" "1528" "1529" "1530"
#> [1531] "1531" "1532" "1533" "1534" "1535" "1536" "1537" "1538" "1539" "1540"
#> [1541] "1541" "1542" "1543" "1544" "1545" "1546" "1547" "1548" "1549" "1550"
#> [1551] "1551" "1552" "1553" "1554" "1555" "1556" "1557" "1558" "1559" "1560"
#> [1561] "1561" "1562" "1563" "1564" "1565" "1566" "1567" "1568" "1569" "1570"
#> [1571] "1571" "1572" "1573" "1574" "1575" "1576" "1577" "1578" "1579" "1580"
#> [1581] "1581" "1582" "1583" "1584" "1585" "1586" "1587" "1588" "1589" "1590"
#> [1591] "1591" "1592" "1593" "1594" "1595" "1596" "1597" "1598" "1599" "1600"
#> [1601] "1601" "1602" "1603" "1604" "1605" "1606" "1607" "1608" "1609" "1610"
#> [1611] "1611" "1612" "1613" "1614" "1615" "1616" "1617" "1618" "1619" "1620"
#> [1621] "1621" "1622" "1623" "1624" "1625" "1626" "1627" "1628" "1629" "1630"
#> [1631] "1631" "1632" "1633" "1634" "1635" "1636" "1637" "1638" "1639" "1640"
#> [1641] "1641" "1642" "1643" "1644" "1645" "1646" "1647" "1648" "1649" "1650"
#> [1651] "1651" "1652" "1653" "1654" "1655" "1656" "1657" "1658" "1659" "1660"
#> [1661] "1661" "1662" "1663" "1664" "1665" "1666" "1667" "1668" "1669" "1670"
#> [1671] "1671" "1672" "1673" "1674" "1675" "1676" "1677" "1678" "1679" "1680"
#> [1681] "1681" "1682" "1683" "1684" "1685" "1686" "1687" "1688" "1689" "1690"
#> [1691] "1691" "1692" "1693" "1694" "1695" "1696" "1697" "1698" "1699" "1700"
#> [1701] "1701" "1702" "1703" "1704" "1705" "1706" "1707" "1708" "1709" "1710"
#> [1711] "1711" "1712" "1713" "1714" "1715" "1716" "1717" "1718" "1719" "1720"
#> [1721] "1721" "1722" "1723" "1724" "1725" "1726" "1727" "1728" "1729" "1730"
#> [1731] "1731" "1732" "1733" "1734" "1735" "1736" "1737" "1738" "1739" "1740"
#> [1741] "1741" "1742" "1743" "1744" "1745" "1746" "1747" "1748" "1749" "1750"
#> [1751] "1751" "1752" "1753" "1754" "1755" "1756" "1757" "1758" "1759" "1760"
#> [1761] "1761" "1762" "1763" "1764" "1765" "1766" "1767" "1768" "1769" "1770"
#> [1771] "1771" "1772" "1773" "1774" "1775" "1776" "1777" "1778" "1779" "1780"
#> [1781] "1781" "1782" "1783" "1784" "1785" "1786" "1787" "1788" "1789" "1790"
#> [1791] "1791" "1792" "1793" "1794" "1795" "1796" "1797" "1798" "1799" "1800"
#> [1801] "1801" "1802" "1803" "1804" "1805" "1806" "1807" "1808" "1809" "1810"
#> [1811] "1811" "1812" "1813" "1814" "1815" "1816" "1817" "1818" "1819" "1820"
#> [1821] "1821" "1822" "1823" "1824" "1825" "1826" "1827" "1828" "1829" "1830"
#> [1831] "1831" "1832" "1833" "1834" "1835" "1836" "1837" "1838" "1839" "1840"
#> [1841] "1841" "1842" "1843" "1844" "1845" "1846" "1847" "1848" "1849" "1850"
#> [1851] "1851" "1852" "1853" "1854" "1855" "1856" "1857" "1858" "1859" "1860"
#> [1861] "1861" "1862" "1863" "1864" "1865" "1866" "1867" "1868" "1869" "1870"
#> [1871] "1871" "1872" "1873" "1874" "1875" "1876" "1877" "1878" "1879" "1880"
#> [1881] "1881" "1882" "1883" "1884" "1885" "1886" "1887" "1888" "1889" "1890"
#> [1891] "1891" "1892" "1893" "1894" "1895" "1896" "1897" "1898" "1899" "1900"
#> [1901] "1901" "1902" "1903" "1904" "1905" "1906" "1907" "1908" "1909" "1910"
#> [1911] "1911" "1912" "1913" "1914" "1915" "1916" "1917" "1918" "1919" "1920"
#> [1921] "1921" "1922" "1923" "1924" "1925" "1926" "1927" "1928" "1929" "1930"
#> [1931] "1931" "1932" "1933" "1934" "1935" "1936" "1937" "1938" "1939" "1940"
#> [1941] "1941" "1942" "1943" "1944" "1945" "1946" "1947" "1948" "1949" "1950"
#> [1951] "1951" "1952" "1953" "1954" "1955" "1956" "1957" "1958" "1959" "1960"
#> [1961] "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969" "1970"
#> [1971] "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980"
#> [1981] "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990"
#> [1991] "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000"
#>
#> $chain
#> [1] "1" "2" "3" "4"
#>
#> $variable
#> [1] "b_specieslevin_cod" "b_specieslevin_fle" "b_phi_specieslevin_cod"
#> [4] "b_phi_specieslevin_fle" "lp__"
pal_diag <- rev(brewer.pal(n = 4, name = "Dark2"))
mcmc_trace(posterior,
pars = c("b_specieslevin_cod", "b_specieslevin_fle",
"b_phi_specieslevin_cod", "b_phi_specieslevin_fle"),
facet_args = list(ncol = 1, strip.position = "left")) +
theme(text = element_text(size = 12), strip.text = element_text(size = 10),
legend.position = "top") +
scale_color_manual(values = alpha(pal_diag, alpha = 0.6))
#> Scale for 'colour' is already present. Adding another scale for 'colour',
#> which will replace the existing scale.
ggsave("figures/supp/levin_mcmc_trace.png", width = 6.5, height = 6.5, dpi = 600)
# Marginal effects of Levin's species covariate
beta_levin_pred <- m_levin_beta %>%
epred_draws(newdata = tibble(species = c("levin_cod", "levin_fle")),
re_formula = NA) %>%
mutate(species = ifelse(species == "levin_fle", "Flounder", "Cod"))
#> mutate (grouped): changed 16,000 values (100%) of 'species' (0 new NA)
p1 <- ggplot(beta_levin_pred, aes(x = .epred, fill = species)) +
stat_halfeye(.width = c(0.8, 0.95), point_interval = "median_hdi", alpha = 0.5,
adjust = 2, position = position_dodge(width = 0.01)) +
coord_cartesian(ylim = c(0, 0.45)) +
scale_fill_brewer(palette = "Dark2", name = "Species") +
theme_light(base_size = 16) +
theme(legend.position = "top") +
guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
labs(x = "Predicted Levin's diversity index", y = NULL
#, caption = "80% and 95% credible intervals shown in black"
) +
NULL
#ggsave("figures/levin_marginal.png", width = 6.5, height = 6.5, dpi = 600)
# Now plot the difference!
beta_levin_pred_cod <- m_levin_beta %>%
epred_draws(newdata = tibble(species = "levin_cod"),
re_formula = NA) %>%
rename(".epred_cod" = ".epred") %>%
ungroup() %>%
dplyr::select(-species)
#> rename: renamed one variable (.epred_cod)
#> ungroup: no grouping variables
beta_levin_pred_fle <- m_levin_beta %>%
epred_draws(newdata = tibble(species = "levin_fle"),
re_formula = NA) %>%
rename(".epred_fle" = ".epred") %>%
ungroup() %>%
dplyr::select(-species)
#> rename: renamed one variable (.epred_fle)
#> ungroup: no grouping variables
beta_levin_pred_wide <- left_join(beta_levin_pred_cod, beta_levin_pred_fle) %>%
mutate(.epred_diff = .epred_cod - .epred_fle)
#> Joining, by = c(".row", ".chain", ".iteration", ".draw")
#> left_join: added one column (.epred_fle)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 8,000
#> > =======
#> > rows total 8,000
#> mutate: new variable '.epred_diff' (double) with 8,000 unique values and 0% NA
p2 <- ggplot(beta_levin_pred_wide, aes(x = .epred_diff)) +
stat_dotsinterval(quantiles = 100) +
geom_vline(xintercept = 0, linetype = 2, color = "tomato", size = 1) +
theme_light(base_size = 16) +
guides(fill = guide_legend(title.position = "top", title.hjust = 0.5)) +
labs(x = "Predicted difference in Levin's diversity index", y = NULL
, caption = "80% and 95% credible intervals shown in black"
) +
NULL
p1 / p2 #+ plot_annotation(tag_levels = "A")
ggsave("figures/levin_species_combined.png", width = 6.5, height = 6.5, dpi = 600)
# This is just very exploratory...
fle_prey_long3 <- fle_prey %>%
mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, area, size_class) %>%
summarise(fle_stomach_content = sum(value)) %>%
arrange(name, size_class, year, area, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, area, size_class, sep = "_")) %>%
group_by(id) %>%
mutate(fle_stomach_content_tot = sum(fle_stomach_content),
fle_stomach_content_prop = fle_stomach_content / fle_stomach_content_tot) %>%
ungroup()
#> mutate: new variable 'size_class' (factor) with 8 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x33, now 34880x19]
#> group_by: 5 grouping variables (name, year, quarter, area, size_class)
#> summarise: now 1,952 rows and 6 columns, 4 group variables remaining (name, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 122 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'fle_stomach_content_tot' (double) with 121 unique values and 0% NA
#> new variable 'fle_stomach_content_prop' (double) with 676 unique values and 2% NA
#> ungroup: no grouping variables
# Now cod
cod_prey_long3 <- cod_prey %>%
mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>%
pivot_longer(15:30) %>%
group_by(name, year, quarter, area, size_class) %>%
summarise(cod_stomach_content = sum(value)) %>%
arrange(name, size_class, year, area, quarter) %>%
ungroup() %>%
mutate(id = paste(year, quarter, area, size_class, sep = "_")) %>%
group_by(id) %>%
mutate(cod_stomach_content_tot = sum(cod_stomach_content),
cod_stomach_content_prop = cod_stomach_content / cod_stomach_content_tot) %>%
ungroup()
#> mutate: new variable 'size_class' (factor) with 16 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> group_by: 5 grouping variables (name, year, quarter, area, size_class)
#> summarise: now 4,080 rows and 6 columns, 4 group variables remaining (name, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'id' (character) with 255 unique values and 0% NA
#> group_by: one grouping variable (id)
#> mutate (grouped): new variable 'cod_stomach_content_tot' (double) with 242 unique values and 0% NA
#> new variable 'cod_stomach_content_prop' (double) with 1,346 unique values and 1% NA
#> ungroup: no grouping variables
# Size based Levin's index
levin_cod <- cod_prey_long3 %>%
drop_na(name) %>%
drop_na(cod_stomach_content_prop) %>%
group_by(id, year, quarter, area, size_class) %>%
summarise(levin = ((1/(number_of_prey$n-1)) * (((1/sum(cod_stomach_content_prop^2))) - 1))) %>%
ungroup() %>%
mutate(size = as.integer(stringr::str_extract(size_class, "\\d+")))
#> drop_na: no rows removed
#> drop_na: removed 48 rows (1%), 4,032 rows remaining
#> group_by: 5 grouping variables (id, year, quarter, area, size_class)
#> summarise: now 252 rows and 6 columns, 4 group variables remaining (id, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'size' (integer) with 16 unique values and 0% NA
levin_fle <- fle_prey_long3 %>%
drop_na(name) %>%
drop_na(fle_stomach_content_prop) %>%
group_by(id, year, quarter, area, size_class) %>%
summarise(levin = ((1/(number_of_prey$n-1)) * (((1/sum(fle_stomach_content_prop^2))) - 1))) %>%
ungroup() %>%
mutate(size = as.integer(stringr::str_extract(size_class, "\\d+")))
#> drop_na: no rows removed
#> drop_na: removed 32 rows (2%), 1,920 rows remaining
#> group_by: 5 grouping variables (id, year, quarter, area, size_class)
#> summarise: now 120 rows and 6 columns, 4 group variables remaining (id, year, quarter, area)
#> ungroup: no grouping variables
#> mutate: new variable 'size' (integer) with 8 unique values and 0% NA
ggplot(levin_cod, aes(size, size_class)) + geom_point()
# Plot Levin
p1 <- levin_cod %>% filter(size > 5 & size < 60) %>%
ggplot(., aes(size, levin)) +
geom_point() +
stat_smooth() +
facet_wrap(~ area, scales = "free", ncol = 1) +
ggtitle("Cod")
#> filter: removed 46 rows (18%), 206 rows remaining
p2 <- levin_fle %>% filter(size > 5 & size < 40) %>%
ggplot(., aes(size, levin)) +
geom_point() +
stat_smooth() +
facet_wrap(~ area, scales = "free", ncol = 1) +
ggtitle("Fle")
#> filter: removed 3 rows (2%), 117 rows remaining
p1 + p2
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# How can I know this is not simply due to samples size? (more samples, more diversity)
# Plot n prey in individual stomachs vs sample size in the year, quarter, ices_rect and size_class
# If there's a strong positive relationship, larger mean number of preys
saturation_fle_prey <- fle_prey %>%
mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>%
pivot_longer(15:30) %>%
filter(value > 0) %>%
group_by(unique_pred_id, year, quarter, area, size_class) %>%
summarise(n_prey = n()) %>%
ungroup() %>%
group_by(year, quarter, area, size_class) %>%
summarise(mean_n_prey = mean(n_prey),
sample_size_n = n()) %>%
ungroup()
#> mutate: new variable 'size_class' (factor) with 8 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x33, now 34880x19]
#> filter: removed 30,910 rows (89%), 3,970 rows remaining
#> group_by: 5 grouping variables (unique_pred_id, year, quarter, area, size_class)
#> summarise: now 1,785 rows and 6 columns, 4 group variables remaining (unique_pred_id, year, quarter, area)
#> ungroup: no grouping variables
#> group_by: 4 grouping variables (year, quarter, area, size_class)
#> summarise: now 120 rows and 6 columns, 3 group variables remaining (year, quarter, area)
#> ungroup: no grouping variables
ggplot(saturation_fle_prey, aes(sample_size_n, mean_n_prey, color = area)) +
geom_jitter(alpha = 0.8, size = 2) +
stat_smooth(method = "lm", se = FALSE) +
facet_wrap(~ size_class, scales = "free") +
scale_color_viridis(discrete = TRUE)
#> `geom_smooth()` using formula 'y ~ x'
saturation_cod_prey <- cod_prey %>%
mutate(size_class = cut(pred_cm, breaks = c(seq(0, 100, 5)))) %>%
pivot_longer(15:30) %>%
filter(value > 0) %>%
group_by(unique_pred_id, year, quarter, area, size_class) %>%
summarise(n_prey = n()) %>%
ungroup() %>%
group_by(year, quarter, area, size_class) %>%
summarise(mean_n_prey = mean(n_prey),
sample_size_n = n()) %>%
ungroup()
#> mutate: new variable 'size_class' (factor) with 16 unique values and 0% NA
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x33, now 57488x19]
#> filter: removed 51,698 rows (90%), 5,790 rows remaining
#> group_by: 5 grouping variables (unique_pred_id, year, quarter, area, size_class)
#> summarise: now 3,229 rows and 6 columns, 4 group variables remaining (unique_pred_id, year, quarter, area)
#> ungroup: no grouping variables
#> group_by: 4 grouping variables (year, quarter, area, size_class)
#> summarise: now 252 rows and 6 columns, 3 group variables remaining (year, quarter, area)
#> ungroup: no grouping variables
saturation_cod_prey %>%
mutate(size = as.integer(stringr::str_extract(size_class, "\\d+"))) %>%
ggplot(., aes(size, mean_n_prey)) +
geom_point() +
stat_smooth()
#> mutate: new variable 'size' (integer) with 16 unique values and 0% NA
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
ggplot(saturation_cod_prey, aes(sample_size_n, mean_n_prey, color = area)) +
geom_jitter(alpha = 0.8) +
facet_wrap(~ size_class, scales = "free") +
stat_smooth(method = "lm", se = FALSE) +
scale_color_viridis(discrete = TRUE) +
guides(color = FALSE)
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
#> `geom_smooth()` using formula 'y ~ x'